//#include <math.h>
//#include <cutil.h>
#include <Particula.h>
#include <vector_types.h>
#include <math.h>

__constant__ float constGas = 1.5;
__constant__ float coeficienteViscosidade = 0.2;
__constant__ float passoIntegracao = 0.001;

__device__ float CalculaKernelME(float3* VBO, int i, int j, float h, int indiceAtual, int resolucaoGridX, int resolucaoGridY, unsigned int* offSet, unsigned int* count, unsigned int* values, float massa, const int numeroParticulas) {
    float massaEsp = 0.0f;
    float kernelPoly = 0.0f;
    for (int a = -1; a < 2; a++) {
        for (int b = -1; b < 2; b++) {
            
            int iTemp = i + a;
            int jTemp1 = j - 1;
            int jTemp2 = j;
            int jTemp3 = j + 1;
            int idCell1 = -1;
            int idCell2 = -1;
            int idCell3 = -1;
            int kInicial1 = -1;
            int kFinal1 = -1;
            int kInicial2 = -1;
            int kFinal2 = -1;
            int kInicial3 = -1;
            int kFinal3 = -1;
            if ((iTemp > 0) && (iTemp < (resolucaoGridY - 1)) && (jTemp1 > 0) && (jTemp1 < (resolucaoGridX - 1))) {
                idCell1 = (iTemp * resolucaoGridX) + (j - 1);
            }
            if ((iTemp > 0) && (iTemp < (resolucaoGridY - 1)) && (jTemp2 > 0) && (jTemp2 < (resolucaoGridX - 1))) {
                idCell2 = (i * resolucaoGridX) + j;
            }
            if ((iTemp > 0) && (iTemp < (resolucaoGridY - 1)) && (jTemp3 > 0) && (jTemp3 < (resolucaoGridX - 1))) {
                idCell3 = (i * resolucaoGridX) + (j + 1);
            }

            if (idCell1 >= 0) {
                kInicial1 = offSet[idCell1];
                kFinal1 = offSet[idCell1] + count[idCell1];
            }

            for (int c = 0; c < numeroParticulas; c++) {
                if ((kInicial1 >= 0) && (kFinal1 >= 0)) {

                    int k = c + kInicial1;
                    if (k < kFinal1) {
                        int indiceParticulas = values[k];
                        if (indiceParticulas != indiceAtual) {
                            float posXAtual = VBO[indiceAtual].x;
                            float posYAtual = VBO[indiceAtual].y;
                            float posXVizinho = VBO[indiceParticulas].x;
                            float posYVizinho = VBO[indiceParticulas].y;
                            float diferencaPosX = posXAtual - posXVizinho;
                            float diferencaPosX2 = diferencaPosX * diferencaPosX;
                            float diferencaPosY = posYAtual - posYVizinho;
                            float diferencaPosY2 = diferencaPosY * diferencaPosY;
                            float distancia2 = diferencaPosX2 + diferencaPosY2;
                            float distancia = sqrt(distancia2);
                            float temp = distancia * 100;
                            float temp2 = h * 100;
                            //distancia = sqrt((VBO[indiceAtual].x - VBO[indiceParticulas].x) * (VBO[indiceAtual].x - VBO[indiceParticulas].x) + (VBO[indiceAtual].y - VBO[indiceParticulas].y) * (VBO[indiceAtual].y - VBO[indiceParticulas].y));
                            //if (diferencaPosX > 0) {
                                massaEsp = posXAtual;
                                //float h9 = h*h*h*h*h*h*h*h*h;
                                //kernelPoly = 315 / (64 * 3.14 * h9) * (h * h - distancia * distancia)*(h * h - distancia * distancia)*(h * h - distancia * distancia);
                                //kernelPoly = 315 / (64 * 3.14 * pow(h, 9)) * pow((h * h - distancia * distancia), 3);
                                //massaEsp += massa * kernelPoly;
                            //}
                        }
                    }
                }
            }
        }
    }

    return massaEsp;
}

/*__device__ float3 CalculaKernelPressao(Particula* particulas, Particula* particulas2, float3* VBO, int i, int j, float h, int indiceAtual, int resolucaoGridX, unsigned int* offSet, unsigned int* count, unsigned int* values, float massa) {
    float3 forcaPressao;
    forcaPressao.x = 0;
    forcaPressao.y = 0;
    forcaPressao.z = 0.0f;
    int idCell = (i * resolucaoGridX) + j;
    float h2 = h*h;

    int kInicial = offSet[idCell];
    int kFinal = offSet[idCell] + count[idCell];

    for (int k = kInicial; k < kFinal; k++) {
        int indiceParticulas = values[k];
        if (indiceParticulas != indiceAtual) {
            float distancia2 = (VBO[indiceAtual].x - VBO[indiceParticulas].x) * (VBO[indiceAtual].x - VBO[indiceParticulas].x) + (VBO[indiceAtual].y - VBO[indiceParticulas].y) * (VBO[indiceAtual].y - VBO[indiceParticulas].y);
            if (distancia2 < h2) {
                float distancia = sqrt(distancia2);
                float3 gradientePressao;
                float mod;
                mod = sqrt(pow(VBO[indiceAtual].x - VBO[indiceParticulas].x, 2) +
                        pow(VBO[indiceAtual].y - VBO[indiceParticulas].y, 2) +
                        pow(VBO[indiceAtual].z - VBO[indiceParticulas].z, 2));
                gradientePressao.x = (VBO[indiceParticulas].x - VBO[indiceAtual].x) / mod;
                gradientePressao.y = (VBO[indiceParticulas].y - VBO[indiceAtual].y) / mod;
                gradientePressao.z = (VBO[indiceParticulas].z - VBO[indiceAtual].z) / mod;

                float kernelPress = (15 / (3.14 * pow(h, 6))) * (pow((h - distancia), 3));
                forcaPressao.x += massa * ((particulas2[indiceParticulas].pressao + particulas[indiceAtual].pressao) / 2 * particulas[indiceParticulas].massaEspecifica) * gradientePressao.x * kernelPress;
                forcaPressao.y += massa * ((particulas2[indiceParticulas].pressao + particulas[indiceAtual].pressao) / 2 * particulas[indiceParticulas].massaEspecifica) * gradientePressao.y * kernelPress;
            }
        }
    }
    return forcaPressao;
}

__device__ float3 CalculaKernelVisc(Particula* particulas, Particula* particulas2, float3* VBO, int i, int j, float h, int indiceAtual, int resolucaoGridX, unsigned int* offSet, unsigned int* count, unsigned int* values, float massa) {
    float3 forcaVisc;
    forcaVisc.x = 0;
    forcaVisc.y = 0;
    forcaVisc.z = 0;
    int idCell = (i * resolucaoGridX) + j;

    int kInicial = offSet[idCell];
    int kFinal = offSet[idCell] + count[idCell];

    for (int k = kInicial; k < kFinal; k++) {
        int indiceParticulas = values[k];
        if (indiceParticulas != indiceAtual) {
            float distancia = sqrt((VBO[indiceAtual].x - VBO[indiceParticulas].x) * (VBO[indiceAtual].x - VBO[indiceParticulas].x) + (VBO[indiceAtual].y - VBO[indiceParticulas].y) * (VBO[indiceAtual].y - VBO[indiceParticulas].y));
            if (distancia > 0 && distancia < h) {
                float kernelVisc = (45 / (3.14 * pow(h, 6))) * (h - distancia);
                forcaVisc.x += massa * ((particulas[indiceParticulas].velocidade.x + particulas2[indiceAtual].velocidade.x) / 2 * particulas[indiceParticulas].massaEspecifica) * kernelVisc;
                forcaVisc.y += massa * ((particulas[indiceParticulas].velocidade.y + particulas2[indiceAtual].velocidade.y) / 2 * particulas[indiceParticulas].massaEspecifica) * kernelVisc;
            }
        }
    }
    return forcaVisc;
}

__device__ void corrigePosicao(float3* posicao, int posicaoParticulaAtual, float XminGrid, float XmaxGrid, float YminGrid, float YmaxGrid){

    if(posicao[posicaoParticulaAtual].x < XminGrid){
        float a = (posicao[posicaoParticulaAtual].x - XminGrid);
        posicao[posicaoParticulaAtual].x = XminGrid + (((a) < 0) ? -(a) : (a));
    }
    if(posicao[posicaoParticulaAtual].x > XmaxGrid){
        float b = (posicao[posicaoParticulaAtual].x - XmaxGrid);
        posicao[posicaoParticulaAtual].x = XmaxGrid - (((b) < 0) ? -(b) : (b));
    }
    if(posicao[posicaoParticulaAtual].y < YminGrid){
        float c = (posicao[posicaoParticulaAtual].y - YminGrid);
        posicao[posicaoParticulaAtual].y = YminGrid + (((c) < 0) ? -(c) : (c));
    }
    if(posicao[posicaoParticulaAtual].y > YmaxGrid){
        float d = (posicao[posicaoParticulaAtual].y - YmaxGrid);
        posicao[posicaoParticulaAtual].y = YmaxGrid + (((d) < 0) ? -(d) : (d));
    }
}*/

__global__ void FindNeighbors(float3* posicoesParticulasT, float3*  posicoesParticulasT2, Particula* particulas, Particula* particulas2, unsigned int* values, unsigned int* keys,
        unsigned int* count, unsigned int* offSet, int resolucaoGridX, int resolucaoGridY, const int numeroParticulas,
        float deltaX, float deltaY, float massaEspecificaInicial, float massa, float h, float3* forcaPressaoVetor,
        float3* forcaAceleracaoVetor, float3* forcaViscosidadeVetor, float XminGrid, float XmaxGrid, float YminGrid, float YmaxGrid, float3 gravidade) {

    int thread = threadIdx.x;
    int block = blockIdx.x;
    int dim = blockDim.x;
    int indiceValue = dim * block + thread;
    int jAtual = 0;
    int iAtual = 0;
    int jMin = 0;
    int jMax = resolucaoGridX;
    int iMin = 0;
    int iMax = resolucaoGridY;
    //int idCell = 0;
    //int kInicial = 0;
    //int kFinal = 0;
    float3 forcaPressao = {0.0f, 0.0f, 0.0f};
    float3 forcaAceleracao = {0.0f, 0.0f, 0.0f};
    float3 forcaViscosidade = {0.0f, 0.0f, 0.0f};
    //float3 gravidade = {0.0f, -9.8f, 0.0f};

    if (indiceValue < numeroParticulas) {
        int posicaoParticulaAtual = values[indiceValue];
        jAtual = (posicoesParticulasT[posicaoParticulaAtual].x / deltaX);
        iAtual = (posicoesParticulasT[posicaoParticulaAtual].y / deltaY);
        //int idCellAtual = (iAtual * numeroColunasGridCell) + jAtual;

        //VERIFICAO DAS EXTREMIDADES.
        /*if (jAtual <= 0) {
            jMin = 0;
        } else {
            jMin = jAtual - 1;
        }

        if (jAtual >= (resolucaoGridX - 1)) {
            jMax = jAtual;
        } else {
            jMax = jAtual + 1;
        }

        if (iAtual <= 0) {
            iMin = 0;
        } else {
            iMin = iAtual - 1;
        }

        if (iAtual >= (resolucaoGridY - 1)) {
            //iMax = iAtual;
            iMax = (resolucaoGridY - 1);
        } else {
            iMax = iAtual + 1;
        }*/

        particulas2[posicaoParticulaAtual].massaEspecifica = 0;

        //particulas2[posicaoParticulaAtual].massaEspecifica += CalculaKernelME(posicoesParticulasT, iAtual, jAtual, h, posicaoParticulaAtual, resolucaoGridX, offSet, count, values, massa, numeroParticulas);
        //CALCULA AS MASSAS ESPECIFICAS.
        //for (int i = -1; i <= 1; i++) {
            //for (int j = -1; j <= 1; j++) {
                int iTemp = iAtual;// + i;
                int jTemp = jAtual;// + j;
                //VERIFICAO DAS EXTREMIDADES.
                if ((iTemp > 0) && (iTemp < (resolucaoGridY - 1)) && (jTemp > 0) && (jTemp < (resolucaoGridX - 1))) {
                    particulas2[posicaoParticulaAtual].massaEspecifica += CalculaKernelME(posicoesParticulasT, iTemp, jTemp, h, posicaoParticulaAtual, resolucaoGridX, resolucaoGridY, offSet, count, values, massa, numeroParticulas);
                }
            //}
        //}

        //CALCALA A PRESSAO NA PARTICULA.
        //particulas2[posicaoParticulaAtual].pressao = constGas * (particulas2[posicaoParticulaAtual].massaEspecifica - massaEspecificaInicial);

        //CALCULA A FORÇA DE PRESSAO.
        /*for (int i = iMin; i <= iMax; i++) {
            for (int j = jMin; j <= jMax; j++) {
                float3 forcaPressaoTemp = CalculaKernelPressao(particulas, particulas2, posicoesParticulasT, i, j, h, posicaoParticulaAtual, resolucaoGridX, offSet, count, values, massa);
                forcaPressao.x += forcaPressaoTemp.x;
                forcaPressao.y += forcaPressaoTemp.y;
                forcaPressao.z += forcaPressaoTemp.z;
            }
        }
        forcaPressao.x *= -1;
        forcaPressao.y *= -1;
        forcaPressao.z *= -1;
        forcaPressaoVetor[posicaoParticulaAtual].x = forcaPressao.x;
        forcaPressaoVetor[posicaoParticulaAtual].y = forcaPressao.y;
        forcaPressaoVetor[posicaoParticulaAtual].z = forcaPressao.z;

        //CALCULA A FORÇA VISCOSIDADE.
        for (int i = iMin; i <= iMax; i++) {
            for (int j = jMin; j <= jMax; j++) {
                float3 forcaViscTemp = CalculaKernelVisc(particulas, particulas2, posicoesParticulasT, i, j, h, posicaoParticulaAtual, resolucaoGridX, offSet, count, values, massa);
                forcaViscosidade.x += forcaViscTemp.x;
                forcaViscosidade.y += forcaViscTemp.y;
                forcaViscosidade.z += forcaViscTemp.z;
            }
        }
        forcaViscosidade.x *= coeficienteViscosidade;
        forcaViscosidade.y *= coeficienteViscosidade;
        forcaViscosidade.z *= coeficienteViscosidade;
        forcaViscosidadeVetor[posicaoParticulaAtual].x = forcaViscosidade.x;
        forcaViscosidadeVetor[posicaoParticulaAtual].y = forcaViscosidade.y;
        forcaViscosidadeVetor[posicaoParticulaAtual].z = forcaViscosidade.z;

        forcaAceleracao.x = forcaPressao.x + forcaViscosidade.x;
        forcaAceleracao.y = forcaPressao.y + forcaViscosidade.y;
        forcaAceleracao.z = forcaPressao.z + forcaViscosidade.z;
        forcaAceleracaoVetor[posicaoParticulaAtual].x = forcaAceleracao.x;
        forcaAceleracaoVetor[posicaoParticulaAtual].y = forcaAceleracao.y;
        forcaAceleracaoVetor[posicaoParticulaAtual].z = forcaAceleracao.z;

        particulas2[posicaoParticulaAtual].velocidade.x = particulas[posicaoParticulaAtual].velocidade.x + passoIntegracao / 2 * forcaAceleracao.x;
        particulas2[posicaoParticulaAtual].velocidade.y = particulas[posicaoParticulaAtual].velocidade.y + passoIntegracao / 2 * forcaAceleracao.y;
        particulas2[posicaoParticulaAtual].velocidade.z = 0;*/

        posicoesParticulasT2[posicaoParticulaAtual].x += 0.05;//posicoesParticulasT[posicaoParticulaAtual].x + passoIntegracao / 2 * particulas2[posicaoParticulaAtual].velocidade.x;
        //posicoesParticulasT2[posicaoParticulaAtual].y = posicoesParticulasT[posicaoParticulaAtual].y + passoIntegracao / 2 * particulas2[posicaoParticulaAtual].velocidade.y;
        //posicoesParticulasT2[posicaoParticulaAtual].z = 0;

        //corrigePosicao(posicoesParticulasT2, posicaoParticulaAtual, XminGrid, XmaxGrid, YminGrid, YmaxGrid);

    }
}

extern "C" void launch_FindNeighbors(unsigned int blocksPerGrid, unsigned int threadsPerBlock, float3* posicoesParticulasT, float3*  posicoesParticulasT2, Particula* particulas, Particula* particulas2, unsigned int* values, unsigned int* keys,
        unsigned int* count, unsigned int* offSet, int resolucaoGridX, int resolucaoGridY, const int numeroParticulas, float deltaX, float deltaY,
        float massaEspecificaInicial, float massa, float h, float3* forcaPressaoVetorDevice, float3* forcaAceleracaoVetorDevice,
        float3* forcaViscosidadeVetorDevice, float XminGrid, float XmaxGrid, float YminGrid, float YmaxGrid, float3 gravidade) {
    // execute the kernel
    FindNeighbors <<<blocksPerGrid, threadsPerBlock >>>(posicoesParticulasT, posicoesParticulasT2, particulas, particulas2, values, keys, count, offSet, resolucaoGridX, resolucaoGridY, numeroParticulas, deltaX, deltaY, massaEspecificaInicial, massa, h, forcaPressaoVetorDevice, forcaAceleracaoVetorDevice, forcaViscosidadeVetorDevice, XminGrid, XmaxGrid, YminGrid, YmaxGrid, gravidade);
}


